Mon. Not. R. Astron. Soc. 000,[TH5](2009) Printed 2 November 2009 (MN KTeX style file v2.2) 



Wavelet-based Faraday Rotation Measure Synthesis 



P. Frick 1 , D. Sokoloff 2 , R. Stepanov 1 , and R. Beck 3 

1 Institute of Continuous Media Mechanics, Korolyov str. 1, 614013 Perm, Russia 

2 Department of Physics, Moscow University, 119899, Moscow, Russia 

3 Max-Planck-Institut fur Radioastronomie, Aufdem Hiigel 69, 53121 Bonn, Germany 



Accepted 2009 .... Received 2009 in original form 2009 



ABSTRACT 

Faraday Rotation Measure (RM) Synthesis, as a method for analyzing multi-channel observa- 
tions of polarized radio emission to investigate galactic magnetic fields structures, requires the 
definition of complex polarized intensity in the wavelength range — <x> < A 2 < oo. The problem 
is that the measurements at negative A 2 are not possible. We introduce a simple method for 
continuation of the observed complex polarized intensity P(A 2 ) into the domain A 2 < us- 
ing symmetry arguments. The method is suggested in context of magnetic field recognition in 
galactic disks where the magnetic field is supposed to have a maximum in the equatorial plane. 
The method is quite simple when applied to a single Faraday-rotating structure on the line of 
sight. Recognition of several structures on the same line of sight requires a more sophisticated 
technique. We also introduce a wavelet-based algorithm which allows us to consider a set 
of isolated structures in the (<p, A 2 ) plane (where (p is the Faraday depth). The method essen- 
tially improves the possibilities for reconstruction of complicated Faraday structures using the 
capabilities of modern radio telescopes. 

Key words: Methods: polarization - methods: data analysis - galaxies: magnetic fields - RM 
Synthesis - wavelets 



1 INTRODUCTION 

Observations of polarized radio emission are the main sources of 
information on magnetic fields of galaxies. The basic idea of mag- 
netic field analysis fro m polarized r adio emission data originates 
in the classica l paper of iBurnl d 19661) (for a later d evelopment see 
ISokoloffetaT] i ll 9981) ). In particular, iBurnl l l 1966b noted that the 
complex polarized intensity P obtained from a radio source is re- 
lated to the Faraday dispersion function F(<p) as 



where k = 2<f>, and the Fourier transform is defined as 



f(k)e' kx dk, f{k) 



)e~ ikx dx. 



(4) 



P(A 2 ) : 



£ 



Fi&e 2 ^ d(j> 



(1) 



F{<f>) is the fraction of radiation with the Faraday depth (j> multiplied 
by intrinsic complex polarization and it is an important emission 
characteristic of interest. Here the Faraday depth (f> is defined by 



0(z) = -O.81 



r 



(2) 



where B\\ is the line-of-sight magnetic field component measured 
in yuG, n e is the thermal electron density measured in crrT 3 and the 
integral is taken from the observer at z = over the region which 
contains both, magnetic fields and free electrons, and z is measured 
in parsecs. Following Eq. (fTJl P is the inverse Fourier transform 
of F. Correspondingly, the Faraday dispersion function F is the 
Fourier transform of the complex polarized intensity: 



F(4>) = -P(k), 

7T 



(3) 



Implementation of multichannel spectro-polarimetry on mod- 
ern radio telescopes provided ob servations of P over a wide range 
of A (e.g. iHaverkorn etai]|2000t) which made the use of Eq. l[3} 
possible. This is the idea of Faraday Rotation Measure Synthesis 
(RM Synthesis) dBrentiens & de Bruvnll2005l) which opened new 
perspectives in inve stigations of magnetic field of galaxies and 
clusters of galaxies {H averkorn et al. 2003 ; Ide Bruvn & Brentiensl 
l2005t lBeckl2009l : Itleald et alj|2009h . 

A key problem of RM Synthesis application is that P is de- 
fined only for A 2 > and in practice can be observed only in a 
finite spectral band. Moreover, the maximum of P in practice can 
be located outside the available spectral band (see e.g. Fig. [Tp). 
Development of robust methods for the reconstruction of F from P 
in a given spectral range becomes crucial for the practical imple- 
mentation of RM Synthesis. 

Fig-EDsho ws results of RM Synthesis app lied to a standard test 
as exploited bv lBrentiens & de Bruvnl f2005l) . Panel (a) shows the 
function F, which includes three real-valued box-like structures, 
panel (b) - the corresponding polarized intensity P (the dashed hor- 
izontal line shows the spectral window 0.6 < A < 0.78 m). We used 
a channel spacing of 8 A = 0.4 cm. Hereafter, F and P are numeri- 
cally evaluated in arbitrary but mutually consistent units. Note that 
F is in general a complex-valued function. Its modulus defines the 
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Figur e 1. RM Synthesis r econstruction of a standard example from 
iBrentiens & de Bruvnl 120051) : (a) initial F(0) which is chosen purely real; 
(b) amplitude of P1(A 2 ); (c) F(<p) reconstructed with whole domain A 2 > 0: 
real part - thin solid, imaginary part - dashed, amplitude - thick solid; (d) 
F(<f>) reconstructed from the data of spectral band 0.6 < A < 0.78 m. The 
spectral window of observations is indicated in panel (b) by horizontal 
dashed line. 



emission and its phase defines the intrinsic position angle. Panel 
(c) shows the result of the straightforward application of the RM 
Synthesis algorithm to the physical range A 2 > 0, while P(A 2 ) is 
set to zero for all negative A 2 . We see that the real part of the re- 
constructed signal is the same as the initial one (except that it has a 
twice lower amplitude), however, the reconstructed signal obtains a 
substantial imaginary part with a shape which is quite remote from 
the real part. This leads to a change of the emission distribution 
and a loss of any information concerning the position angle (apart 
from the central point of the emission region, where the position 
angle correctly is zero). In the context of ch aotic magnetic fields in 
galax y clusters this loss is less important dde Bruvn & Brentienj 
1 20051) , but in galactic magnetic field studies it becomes crucial 
because the intrinsic position angle determines the orientation of 
the regular magnetic field component perpendicular to the line of 
sight. Fig. [TJi shows that the reconstruction becomes much more 
difficult if we restrict the data to a relatively narrow spectral band 



0.6 < A < 0.78 m. We see that even the sign of the reconstructed 
real part can be wr ong. In that case the algorithm for finite spectral 
band introduced bv lBrentiens & de Bruvnl d2005l) was used. 

A general message obtained from Fig.[T]is that in order to en- 
visage possible ways to get a practical implementation of RM Syn- 
thesis one has to include some additional information based on the 
nature of the physical phenomena which provide the Faraday ro- 
tation. Here we concentrate our efforts on the problems associated 
with missing P(A 2 ) for A 2 < 0. 



2 IMPROVING THE RM SYNTHESIS ALGORITHM 

The complex-valued intensity of polarized radio emission for a 
given wavelength 



P(A 2 ) : 



r 



e{z)e 2i x {z) e 2i ^ )A2 dz, 



(5) 



is defined by the emissivity e and the intrinsic position angle x 
along the line of sight. Here z is the distance from observer to a 
point in the emitting region; the integral is taken over the whole 
emitting region. If the Faraday depth is a monotonic function of z 
(which means that z is a single- valued function of (p), we can define 
the Faraday dispersion function as a function of Faraday depth 



(6) 



In the ideal case, reconstructing the Faraday dispersion function F 
from l[3} and knowing the Faraday depth for any z, one can de- 
rive the characteristics of radio emission (e and^-) along the line of 
sight. They can be used as a tomography in order to derive some 
characteristics of the magnetic field distribution from F. The task of 
RM Synthesis is much more modest and concerns the reconstruc- 
tion of the Faraday function from the observed polarized emission 
which itself is already a complicated problem. 

Let us consider a physically motivated simple example, i.e. P 
produced by a two-layer system , to isolate and overcome the short- 
comings of the RM Synthesis technique. Each layer contains a ho- 
mogeneous magnetic field which has non-vanishing line-of-sight 
and perpendicular components. Both layers are thought to be emit- 
ting and rotating polarized radio waves. The corresponding F(if>) 
is shown in Fig. 2a. It is important for the discussion below that 
the analyzed signal has non-vanishing real and imaginary parts. 
The absolute value of F((f>) indicates how much polarized emission 
comes from a region with Faraday depth (f> and its phase gives the 
intrinsic position angle (about 13° and 31°) of the emission. Just to 
illustrate the variety of possible situations, we choose two different 
shapes of the slabs, i.e. one slab with sharp boundaries and one with 
a Gaussian shape. 

The result of the straightforward application of RM Synthesis 
where the integral is taken over the physically admissible region 
A 2 > is shown in Fig.[2p. RM Synthesis reproduces to some extent 
the absolute value of the signal, but fails to reproduce its phase. 
A naive interpretation of this result could be that field reversals 
occur in each layer, but is obviously incorrect. In the same figure 
we show the result of F((f>) reconstruction within the spectral band 
0.6 < A < 0.78 m (panel c). Then both structures become diffuse 
with a more or less arbitrary phase. The last panel illustrates what 
happens if the upper wavelength boundary will be extended up to 
A = 2.5m (as expected for the Low Frequency Array (LOFAR) 
and the Square Kilometre Array (SKA) telescopes). This extension 
essentially improves the recognition of the sharp structure (the right 
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Figure 2. Standard RM Synthesis for a test Faraday dispersion function, (a) 
Original test function which includes one Gaussian and one box structure. 
Reconstructions: (b) using the whole domain A 2 > 0, (c) using the window 
0.6 < A < 0.78 m, (d) using the window 0.6 < A < 2.5 m. Real part - thin 
solid, imaginary part - dashed, amplitude - thick solid. 



Figure 3. RM Synthesis for the test from Fig|2]using the extension of P(A 2 ) 
in the domain A 2 < defined by (7). The parameter <f>o is adjusted to the 
position of the left structure (a,b) or right structure (c,d). The whole domain 
of A is used in panels (a,c) and the spectral window 0.6 < A < 0.78 m in 
panels (b,d). Real part - thin solid, imaginary part - dashed, amplitude - 
thick solid. 



one in the figure) but almost does not affect the reconstruction of 
the left (Gaussian) structure. 

To avoid the non-uniqueness in the Faraday dispersion func- 
tion reconstruction, some additional information (or hypothesis) is 
required. We suggest to improve the above reconstruction by some 
constraint concerning the possible symmetry of an isolated object. 

Suppose that the expected objects are mainly galactic disks 
with magnetic fields believed to be symmetric with respect to the 
galactic equator. Then the desired F should be even with respect 
to the center of the given object. Therefore, we consider each max- 
imum of the reconstructed F{<p) separately and prescribe that the 
continuation of P(A 2 ) to the region of A 2 < has to be chosen in a 
way which makes F((f>) symmetric with respect to the point (f> = (f> , 
where (f>o is the position of the maximum under consideration. This 
means that F(2(f> -(/)) = F((f>) and using the shift theorem one gets 

P(-A 2 ) = exp(-4i<f> Q A 2 )P(A 2 ). (7) 



The antisymmetric case can be considered as well with slight 
change in the algorithm: Eq. {7]l changes to P(-A 2 ) = 
-exp(-4i(f> A 2 )P(A 2 ). 

Fig. [3] shows the results of reconstruction of the same test but 
following the suggested continuation. The test function includes 
two objects, while the algorithm includes only one parameter <po. 
Firstly, we performed the continuation adjusting (f> to the position 
of the left object (panel a). Then the method gives realistic result 
for this object. The reconstructed structure has no apparent inter- 
nal field reversal and the ratio of real and imaginary parts of F(0), 
i.e. the phase, is correctly reproduced. Position angles are restored 
with the accuracy of 3°. Of course, the result for the other layer, 
i.e. the second maximum of \F((f)\ in Fig. [3] remains false. Panel 
(b) shows what happens if the range of A covered by the observa- 
tion is reduced to 0.6 < A < 0.78 m. Instead of one peak one gets 
a sequence of peaks, which is a usual result for a Fourier recon- 
struction using a narrow spectral window. The suggested procedure 
does not suppress the sidelobes in the standard Rotation Measure 
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Figure 4. Wavelet-based RM Synthesis for the test from Fig|2] The modulus of wavelet coefficients on the (a, b) plane (panels a,c,e) and the result of 
reconstruction (panels b,d,f) for whole domain of A (panels a,b) and the windows 0.6 < A < 0.78 m (panels c,d) and 0.6 < A < 2.5 m (panels e,f). Real part - 
thin solid, imaginary part - dashed, amplitude - thick solid. 



Spread Function (RMSF) dHeald et al.l2009» but corrects the phase 
within the main central peak. Of course, the amplitude of each peak 
is much less than the amplitude of the peak in panel (a), however, 
the ratio of real and imaginary parts of F((p) in the central peak 
remains realistic. If the parameter (p is chosen following the posi- 
tion of the second object the method gives a correct reconstruction 
for the right layer and fails to reproduce the left one. An obvious 
shortcoming of the method exploited is its local nature: We obtain a 
realistic shape of a chosen maximum and ignore what happens with 
the other one. A natural extension is to apply the recommendation 
of Eq. locally to each maximum. This extension brings the idea 
of wavelets into consideration. 



3 RM SYNTHESIS AND WAVELETS 

Wavelet transform presents a kind of "local" Fourier transform, 
allowing us to isolate a given structure in physical space and the 
Fourier space. Let us define the wavelet transform of the Faraday 
dispersion function F(<p) as 



w F (a, b) 



\a\ I 



dip, 



(8) 



where if/(if>) is the analyzing wavelet, a defines the scale and b de- 
fines the position of the wavelet. Then the coefficient w F gives the 
contribution of corresponding structure into the function F. 

The func tion F can be recon structed using the inverse trans- 
form (see, e.g. lDaubechiesI i ll 9921) ) 



1 f ? ld>-b\ dadb 
The reconstruction formula exists under condition that 



\k\ 



dk < 



(10) 



Here iff(k) = J tj/{(p)e ,k ^d<f> is the Fourier transform of the analyzing 
wavelet if/((f>). 



Let us emphasize that the inverse formula is usually writ- 
ten for real signals. Then the scale parameter a is positively defined 
and the integral is taken for < a < oo. In the case of a complex- 
valued function, the range of a can be limited by positive values 
a > by taking a real analyzing wavelet if/(x). In general case of a 
complex- valued function and a complex wavelet, the scale param- 
eter a should be extended into the domain of negative values (like 
wave numbers in Fourier space). 

For the sake of definiteness, we use as the analyzing wavelet 
the so-called Mexican hat if/(if>) = (1 - (j> 2 ) exp(-0 2 /2). The 
wavelet is real, however, the function P is complex, so that the 
wavelet coefficients w F are complex as well. For the chosen wavelet 
Wf(-a,b) = w F (a, b) and = 1. 

Using the definition of the wavelet transform (0 and relation 
we can directly define the wavelet decomposition of the Faraday 
dispersion function from the polarized intensity P(A 2 ) 



w F (a,b) = i J P(A 2 )e- 2ibAl <(/' (-2aA 2 ) 



dA 1 



(11) 



Note that in the case of real F the problem of negative A 2 can 
be solved using progressive wavelets, whose Fourier image is lo- 
calized in the domain of positive wave numbers. Thus using this 
kind of wavelets one avoids the problem of the P(A 2 ) continuation 
in the domain A 2 < 0. 

For the general case, we divide Eq. I ll 11 1 in two parts w F (a, b) = 
W-(a, b) + w+(a, b), where 



w_(a, b) 



w+(a, b) 



u 

i J P(A 2 )e- 2,M2 ^ (-2aA 2 )dA 2 , 
~ J P(A 2 )e- 2,bx2 {p' {-2aA 2 )dA 2 . 



(12) 



(13) 



We propose the following algorithm: Firstly, knowing P(A 2 ) 
for A 2 > we calculate the coefficients w+(a,b) and we recognize 
the dominating structures in the map |w+(a, b)\. The coordinate b of 
the corresponding maximum gives us the value of where upper 
index i indicates the number of the structure. Then we reconstruct 
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Figure 5. The intrinsic position angle x f° r the test from Fig[2ja) has been 
reconstructed with standard (dashed) and wavelet-based (thin solid) RM 
Synthesis. The whole domain of A is used. Thick solid lines show initial x 
in the location of the both structures. 

the coefficients vv_(a, b) following the idea of Eq. Q, but refor- 
mulated for the local domain in wavelet space (a, b). Namely, we 
define 

w_(a, b) = w + (ci, 2f {a, b) - ft) , (14) 

where the parameter (f>' (a,b) for the given point (a, b) is chosen 
according to the structure i which dominates in its vicinity. 

Now we apply the suggested algorithm to the test function 
from Fig. [2] The map b)\ presented in Fig. [4^ demonstrates 

two well-defined structures. The ^-coordinates of the maxima are 
taken as (f>' Q . The result of the reconstruction (see Fig. |4p) shows 
that the method reproduces the amplitude and phase of F(0) for 
both layers. The reconstruction here is performed using P(A 2 ) for 
the whole range A 2 > 0. The comparison of the reconstructed posi- 
tion angle using standard and wavelet-base RM Synthesis is shown 
in Fig. [5] The suggested algorithm gives correct value for;f within 
both emission regions. Panels (c,d) show what happens for the re- 
construction using the spectral window 0.6 < A < 0.78 m. One 
can see the wavelet map is empty in its substantial part a > 2, 
however, the structures remain well-recognizable (panel c). The re- 
constructed F contains several oscillations in domains related to 
both layers. The amplitude of each oscillation becomes much lower 
than that in panel (b), however, the ratio of the real and imagi- 
nary parts in the central maxima remain correct. The third cou- 
ple of panels shows the reconstruction within the extended window 
0.6 < A < 2.5 m. This extension allows one to keep the horn-like 
structures in the bottom of the wavelet plane (panel e) which pro- 
vide the reconstruction of sharp boundaries of the box-like structure 
(panel f). 
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netic field strength is supposed to have a maximum in the equatorial 
plane. 

The suggested method is quite simple when applied to a sin- 
gle structure on the line of sight. Recognition of several structures 
on the same line of sight requires a more sophisticated technique. 
The problem of structure separation is resolved using the wavelet 
decomposition. A simple test example demonstrates the applicabil- 
ity of this method. The polarization angle reconstruction is sig- 
nificantly improved over the standard technique. The wavelets can 
be useful to also overcome some other problems of RM Synthe- 
sis, related to the multi-band structur e of the observational do main 
in /l-space, noise filtration, etc (e.g. iFrick et alii 19971 . 120011) . The 
method essentially improves the possibilities for reconstruction of 
complicated Faraday structures using the capabilities of modern ra- 
dio telescopes. 

Finally note that our simple examples illustrate that the ex- 
tension of the observational band into the long-wavelength domain 
is helpful for the recognition of structures with sharp boundaries, 
while the short-wavelength domain is crucial for the reconstruction 
of smooth structures. 
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4 CONCLUSIONS 

The development of multi-channel observations of polarized ra- 
dio emission opens promising perspectives in the understanding 
of cosmic magnetic fields on galactic and intergalactic scales. The 
first fruitful applications of RM Synthesis suggested in this con- 
text include the recogn ition of local struct ures in the Milky Way 
dHaverkorn et al.l2003h , clu sters of galaxies dde Bruyn & Brentjens! 
120051) and spiral galaxies l lHeald et"al] 120091) . However, in gen- 
eral the RM Synthesis algorithm contains a fundamental prob- 
lem emerging from the fact that the reconstruction formula re- 
quires the definition of complex polarized intensity in the range 
-oo < A 2 < oo. In this paper we introduce a simple method for con- 
tinuation of observed complex polarized intensity P(A 2 ) into the 
domain of negative A 2 < 0. The method is suggested in context 
of magnetic field recognition in galactic disks, for which the mag- 



